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Summary The new multi-grid (or adaptive) pseudospectral element method has been carried out 
for the solution of incompressible flow in terms of primitive variable formulation. The desired features 
of the proposed method include (1) the ability to treat complex geometry; (2) the high resolution 
adapted in the interesting areas; (3) the minimal working space; and (4) effective under the multiple 
processors working environment. 

The approach for flow problems, complex geometry or not, is to first divide the computational 
domain into a number of fine-grid and coarse-grid subdomains with the inter-overlapping area. Next, 
implement the Schwarz alternating procedure (SAP) to exchange the data among subdomains, where 
the coarse-grid correction is used to remove the high frequency error that occurs when the data 
interpolation from the fine-grid subdomain to the coarse-grid subdomain is conducted. The strategy 
behind the coarse-grid correction is to adopt the operator of the divergence of the velocity field, which 
intrinsically links the pressure equation, into this process. The solution of each subdomain can be 
efficiently solved by the direct (or iterative) eigenfunction expansion technique with the least storage 
requirement, i.e., O(N^) in 3-D and 0(iV 2 ) in 2-D. 

Numerical results of both driven cavity and jet flow will be presented in the paper to account for 
the versatility of the proposed method. 
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1 Introduction 


Due to the advance of numerical techniques, numerous CFD algorithms have been developed to 
pursue the hard-to-approach flow problems. Nevertheless, numerical algorithms should have desired 
features of (1) the ability to deal with the variety of geometrical shapes; (2) arbitrary layout of dense 
grid points in the interesting areas; (3) the minimal working space; and (4) the low computational 
time to achieve such a goal. The development of a pseudospectral element method in these areas is 
our major concern. 

One of the improvements in the area of feature (2) is the multi-grid technique, which has long been 
advocated by the finite-difference method [1, 2]. On the same computational domain, a sequence of 
uniform grids are employed to accelerate the convergence of iterative methods. The work rests on the 
“standard coarsening,” i.e., doubling the mesh in each direction from one grid to the next coarsest 
grid and also smoothing the residual to the next coarse grid (restriction). Solve the problem on the 
coarse grid (low frequency domain) and the coarse-grid correction transfers back (prolongation) to 
the fine grid (high frequency domain) to gain rapid convergence. The technique developed so far, 
even with the inclusion of an adaptive scheme, is still limited to the simple complex geometry with 
uniform grids in the Cartesian coordinates, but is less for the non-uniform grids in the curvilinear 
coordinates. - ^ 

The SAP iterative scheme has been successfully applied by the pseudospectral element method 
to those (simple complex) configurations where the overlapped grids are located at the same places 
[3, 4]. Here we refer to such cases as a single-grid SAP method because no error is involved during the 
data interpolation process. But under some circumstances, due to the complexity of the geometrical 
configuration such as possible layout of mixed types of grids (Cartesian, “0” or “C”) or the necessity 
of applying adaptive fine grids for high resolution in one area and coarse grids for less resolution 
in others, the overlapped grids cannot be collapsed at the same position. Careful treatment on the 
overlapped grids by the SAP iterative scheme to eliminate the high frequency error due to the data 
interpolation will be the main objective in this paper. On the other hand, the question arises of 


fiow the continuity equation is satisfied in the overlapping area (including the interfaces between the 
fine-grid and coarse-grid subdomains) when solving the incompressible Navier- Stokes equations in 


primitive variable form. It reflects the fact that the boundary conditions for the pressure should link 
the incompressibility constraint in some respects. Extension of single-grid SAP to the multi-grid SAP 
domain decomposition method to overcome the above-mentioned difficulties will also be addressed. 

The paper consists of five additional sections. Section 2 derives a primitive variable form of the 
Navier-Stokes equations. Section 3 discusses the multi-grid SAP domain decomposition method. 
Section 4 presents numerical results of proposed 2-D problems, and the final section provides the 
conclusions. 


2 Primitive Variable Formulation 


In tensor notation, the time-dependent Navier-Stokes equations in dimensionless form can be de- 
scribed as 


dui dui dp 1 d 2 Ui 

dt + U * dxj dxi Re dxj 


dui 

dxi 


- 0 


(la) 

(lb) 
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Here u* is the velocity component and Re is the Reynolds number. 

The method applied to solve the Navier-Stokes equations is Chorin’s [5] splitting technique. Ac- 
cording to this scheme, the equations of motion read 


dui dp 
dt ' dxi 


= F t 


( 2 ) 


where Fi = —Uj du{/ dxj+ 1 /Re d 2 U{/ dx 2 -. 

The first step is to split the velocity into a sum of predicted and corrected values. The predicted 
velocity is determined by time integration of the momentum equations without the pressure term 

u " +1 = < + A tF? (3) 


The second step is to develop the pressure and corrected velocity fields that satisfy the continuity 
equation by using the relationships 


* . dp 
A( _ 


,"+I — ,-7" +1 


d < +1 

dxi 


= 0 


(4a) 

(4b) 


Here the superscript n denotes the n-th time step. Note that the size of a stable time-step At can be 
increased by using an adaptation of Runge-Kutta techniques [6] for the high Reynolds number and 
the Stokes solution for the low Reynolds number [3], respectively. 

An equation for the pressure can be obtained by taking the divergence of Eq. (4a). In view of Eq. 
(4b), it governs 

d 2 P _ /'c'l 

dx\ At dxi 

Note that whenever solving Eq. (5) the identity of Eq. (4a) should be utilized to absorb the given 
boundary conditions of the velocity components [7]. 

If p satisfies Eq. (5), then u n+1 does indeed satisfy Eq. (4b). The solution of the pressure equation, 
Eq. (5), is the most computationally expensive step, while in Cartesian coordinates it can be directly 
solved numerically by the separation of variables [7]. Eq. (5) is of the general form, 


Lp= S 


( 6 ) 


for some linear operator L on some finite dimensional vector space. The properties of the operator 
L depend on the methods chosen to represent the fields and their derivatives. 

Let the pressure p and source term S in Eq. (6) be expanded in a series of eigenfunctions such 
that 


p = EXpEY T EZ T (7a) 

S = EXSEY t EZ t . (7b) 

then the solution of three-dimensional pressure Eq. (5) can be reduced to the simplest algebraic form 

(oj + fij + 7 k)pi,j,k — Sij t k (8) 
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where a*, /3j and 7 * are the eigenvalues with respect to the discrete derivative matrices of the linear 
operator, L , and EX, EY, EZ are the corresponding eigenvectors associated with each eigenvalue. 
However, eigenvalues may not be real due to the complexity of an operator L. Without putting any 
restriction on eigenvalues, complex eigenvalues and their associated eigenvectors are permitted if the 
pressure gradient at the imaginary part vanishes. This is true because only the pressure gradient 
drives flow instead of the pressure itself. However, the effort of the matrix multiplication will be 
increased by a factor of four if all the calculations of Eqs. (7) are performed by the purely complex 
variables. Fortunately, not exceeding a factor of two will be reached if one takes advantage of (1) the 
purely real part of eigenfunctions for matrix multiplication; (2) the source term S being real; and (3) 
choosing only the real part of pressure as the pressure solution. The way for (1) includes reordering 
the eigenfunction into two parts: real versus complex, and similarly for the eigenvalues. 

The iterative preconditioned method for the solution of pressure in the curvilinear coordinate 
system can be found in [8]. Note that if there are N degrees of freedom in each direction the overall 
memory required for finding the solution to the pressure equation in three dimensions is 0(N '*). This 
is the same type of maximal storage efficient scaling that we have for the velocity field. 

Viewing the solution of the Navier-Stokes equations by the splitting method, two steps account 
for most of the run time, predicted velocity and the pressure solution. The bulk of these two steps 
can be concisely described in terms of dot products and matrix multiplication between subsets of 
array. Importantly, no data dependency occurs when running programs on parallel machines. 


3 Domain Decomposition with Multi-Grid SAP 

The solution of the Navier-Stokes equations via the domain decomposition approach consists of first 
dividing the computational domain into a number of subdomains with inter-overlapping areas, where 
the grids inside the overlapping area may or may not be located at the same places. Next imple- 
ment the SAP for exchanging data among subdomains, i.e., solving the problem on each subdomain 
separately and then updating the velocity field on the overlapped interfaces. The advantages of this 
approach include (i) less memory access, local rather than global, and (ii) easy treatment of complex 
geometry. 

In addition to the Lagrangian constraint between the pressure and velocity fields, the noncoinci- 
dent overlapped grids in the inter-overlapping areas among subdomains even enhance the difficulty 
of applying the multi-grid technique. However, the idea of “coarse-grid correction” is still effective 
to reduce the high frequency error from the interpolated residual of the fine-grid subdomain. The 
strategy behind the coarse-grid correction process is to adopt the idea proposed by Thompson and 
Ferziger [9] and is modified as 

V c • Uc - V c • (#.,) = j/( r, - V, ■ u,) (9) 

Here V c - represents the operator of divergence on the coarse-grid subdomain. 1/ is an interpolation 
operator from the fine-grid subdomain / to coarse grid subdomain c, and u is the velocity component. 
rj is simply the result of the divergence of the velocity field which should be set to zero. The left hand 
side of Eq. (9) is the difference between the coarse-grid operator acting on the coarse-grid subdomain 
and the coarse-grid operator acting on the interpolated fine-grid subdomain (which is held fixed). 
The right hand side of Eq. (9) is the interpolated residual of the fine-grid subdomain. It is obvious 
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that once the solution of the fine-grid subdomain has been found the residual will be zero (exactly 
satisfying the pressure equation), and it also implies 

u c = l/u/ (10) 

When the residual is non-zero, Eq. (9) acts as a forcing term for the coarse-grid correction to transfer 
the correction of the velocity field back to the fine-grid subdomain, i.e., 

u n ™ = u °j d + r f {u c - 1 l u° ld ) (ii) 

This is vital for the success of the scheme. Changes in the velocity field are transferred back to the 
fine-grid subdomain rather than the velocity field itself. Notice that when the overlapped grids in 
the overlapping areas are collapsed at the same places the interpolation operator // automatically 
becomes a unitary matrix. 

The multi-grid SAP iterative solution of the incompressible Navier-Stokes equation in primitive 
variable form for a driven cavity flow sketched in Fig. 1 is summarized by the following algorithm. 

1. First assume u n+1 on AB. Usually u n will be a good initial guess. 

2. Solve fine-grid domain II employing the boundary conditions derived from the divergence of 
the velocity field, including on AB, where the pressure solution is directly obtained by the 
eigenfunction expansion method. 

3. With the interpolated solution of u n+ * from step (2) on domain IIIC I, solve coarse-grid domain 
I employing the same type boundary conditions including on CD to update u n+1 on domain 
III C II by the coarse-grid correction process. 

4. Repeat steps (2) & (3) until the velocity u n+1 on AB, CD does not change. 

In order to guarantee that consistent values of velocity (or pressure gradient) be generated in the 
overlapping domains III, satisfying Eq. (10), the divergence of the velocity field V • u needs to be 
actually computed in whichever domain I or II is counted [4]. Since u on domains III is not known 
a priori, the divergence of the velocity field is only set to zero at the first SAP iteration for step (2). 
According to this approach, the continuity equation is satisfied on domains II (including III C II) 
and I (excluding III C I), which is revealed from Eq. (9) that the continuity equation is only satisfied 
on the fine-grid domain II. More specifically, the issue of how to satisfy the continuity equation along 
the interfaces of fine-coarse grid domain can be easily resolve d by the proposed approach, namely, 
V • u on AB satisfied on the coarse-grid domain I, and V • u on CD satisfied on the fine-grid domain 
II. However, the error index of the continuity equation on domain III C I will indicate how good the 
interpolation is (affected by the layout of overlapped grids) and whether any steep change of flow 
field exists. 

Three main issues occurring in the overlapping area between the fine-grid and coarse-grid subdo- 
mains one might often encounter are how to (1) efficiently implement the interpolation; (2) adequately 
represent the predicted velocity; and (3) explicitly impose the global mass conservation. Each will 
be addressed separately. 
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3.1 Data interpolation 

Finding the image (£,77), -1 < £ < 1,-1 < 17 < 1, of a collocation point (x ,y) from the fine-grid 
subdomain II mapped onto the coarse-grid subdomain I (or vice versa) is first determined by using 
the two-dimensional Lagrange interpolation to seek its corresponding position falling into an element 
on the coarse-grid subdomain that contains (M + 1) x (N -f 1) collocation points, = cosm/M(i — 
0, M), r]j = cosnj /N(j = 0, N ), such that 

M N 

771=0 n=0 
M N 

y= E E K«T m (t)T^) 

771=0 71=0 

where T m denotes the mth order Chebyshev polynomials. Unknown expansion coefficients, a mn , b mn , 
can be easily obtained by the prescribed points ( x , y) on the coarse-grid subdomain I through 

M N 

x {£u r Jj) = ^ a mnT m (£i)T n (r)j) (13a) 

m=0 n=0 
MN 

y{£hVj) ~ bmnTm(£i)T n ('ni) (13b) 

771=0 n=0 

With a given point ( x,y ) in the physical space of fine-grid subdomain II, its image (£,?/) on the 
coarse-grid subdomain I can be iteratively solved by the Newton-Raphson method. Once the one- 
to-one correspondence between the fine-grid and coarse-grid subdomains has been established, the 
equation required to generate a function y) on the fine-grid subdomain interpolated from the 
coarse-grid subdomain, now becomes 

M N 

4‘ , (x, y ) = ee m) (M) 

i=0 j= 0 

where <f> c {£i,rij) denotes the function value at the collocation point rjj) on the coarse-grid subdo- 
main, and lVj(£), Nj(r) ) are the shape functions defining the geometry of the element on the coarse-grid 
subdomain, whose expressions are 


(12a) 

(12b) 


M 

m)= Er-nffXUf,-) (is*) 

m= 0 
N 

l) = E r n(»7)T„(77i) (15b) 

71=0 

where the matrices T m (£) and f m (£) are the Fourier cosine series and their inverse [7]. Note that the 
shape functions N^), Nj(rj) satisfy the Kronecker-delta property, i.e., N t (£ m ) = 6 im ,Nj(T ] n ) = S jn . 

Be aware that it requires much less effort to perform the data interpolation if the one-to-one 
correspondence for the shape functions between subdomains can be stored (once and for all). Also 
the cost for such additional memory is negligible compared to that declared by a single variable. 
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3.2 Predicted velocity 

Since the predicted velocity in the overlapping area generated from Eq. (3) by the fine-grid subdomain 
is slightly different from that obtained by the coarse-grid subdomain, how to control the predicted 
velocity in order to keep the error index, £2 norm of u> = ||u c — J/u/|| minimal, is of great importance. 
Numerical experiments suggest that the following dynamic relationship 

u/ = (1 — u>“) u/ + w a I c f u c (16) 

gives the best fit. Here the exponent a is chosen as 0.4 for various tested problems. 


3.3 Global flow rate 

For the inflow-outflow problems the coarse-grid velocity field interpolated from the fine-grid subdo- 
main may not exactly satisfy the global mass conservation, and a slight adjustment to the velocity 
field is imperative. A common- used formula will meet such a requirement, i.e., 


u 


/ ..exact 




•dA 


u i t 

/»*“■ 


dA 


(17) 


4 Results and Discussion 

For the numerical test of the driven cavity flow problem, layout of elements (6 points per element) 
in the fine- grid and coarse-grid subdomains at the Reynolds number of 400 and 100 are displayed 
in Figs. 2a and 2b, respectively. The overlapping area is not explicitly shown in the figures, but 
just imagine the extension of one more element from the coarse-grid subdomain into the fine-grid 
subdomain. The layout of elements is in accordance with the requirement to resolve the steep 
changes inside the boundary layers. When exchanging the data through the interpolation in the 
inter-overlapping area, the high frequency error introduced by the fine-grid subdomain will pollute 
the results throughout the whole computational domain. It can be simply proved by checking the 
error index, £2 norm of “> = K- J/ u /ll> in the overlapping area, w will increase with marching 
in time domain, and eventually become an unreasonably big number under which the solution does 
not exist. With the multi-grid SAP approach, both results produce O(10 -4 ) for u, instead. When 
comparing the streamline plots, rp m i n = -0.1055 for Re = 100 and = -0.1163 for Re = 400, with 
the most accurate results of Ghia [10], good agreements can be observed in Fig. 3. 

For the inflow-outflow jet problem, a nozzle is designed to gain a high speed fluid with a smooth 
change of the convergent channel. The configuration of jet flow is plotted in Fig. 4. A jet emanating 
from the nozzle with an aspect ratio H/D — 144 (the width of tank to nozzle) is used to understand the 
turbulent characters through the direct numerical simulation. With a strong stratification imposed 
in the vertical direction the two-dimensional turbulent flow calculation will be a good approximation 
to the three-dimensional case. The calculation is carried out up to the Kolmogoroff length scale 
where the energy transferred from the large scales is in equilibrium with the energy dissipated in the 
smallest scale by the molecular viscosity. Certainly, for the purpose of direct numerical simulation 
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the Reynolds number should not be large so that the machine can still handle the huge number of 
points required for the resolution of different length scales. The computational domain is decomposed 
into three subdomains: the upstream nozzle where the inflow is developed to gain a high speed, the 
immediate downstream from the exit of the nozzle where the high speed jet is discharged into the 
tank, and the far downstream (fine grids) where a well- developed turbulent flow can be traced. 

Let us first check the error index u> for the inflow-outflow jet problem without using the multi-grid 
SAP technique. The u> around O(10~ :i ) seems all-right at Re = 100 initially, but the onset of noise 
starts to destabilize the downstream flow field at the Reynolds number of 250 and u increases up 
to O(10 -2 ). That clearly demonstrates the high frequency polluting that results on the fine-grid 
subdomain, but the noise can be totally removed by using the multi-grid SAP technique. Fig. 5 
depicts the streamline plot of jet flow at Re = 250. During the time evolution of the jet flow, the 
symmetry of the jet front will not be distorted at the early stage (Fig. 5a) until the phase speed of 
vortex shedding (due to flow instability) travels faster than that of the jet front. As seen in Fig. (5b), 
a pair of vortices adjacent to the jet front persisting throughout the course represent the extrusion 
of the jet into the ambient fluid. Once the jet front is caught up by the incoming travelling waves, 
the energy transferred by the vortex shedding, in terms of the cascade process from the highest at 
the nozzle exit (high shedding frequency) to the lowest at the jet front (low shedding frequency), 
splits into two parts, one for the jet front to push against the ambient viscous resistance, another 
for the vertical motion. The intensity of vertical motion behind the jet front is gradually enhanced 
as visualized by the splitting streamlines, and their patterns move backward toward the exit of the 
nozzle where a distinct pair of vortices exist. The appearance of similar pairs of vortices can also be 
confirmed by the experiment at the high Reynolds number [11]. 


5 Conclusions 

The solution of the Navier-Stokes equations in a primitive variable form has been solved by the pseu- 
dospectral element method via the multi-grid domain decomposition technique. The computational 
domain is divided into a number of simple subdomains with the inter-overlapping zone, of which the 
fine grids (or fine-grid subdomain) are used in the areas with the steep change of flow field while the 
coarse grids (or coarse-grid subdomain) are used in the others. During the data exchange among 
subdomains, the coarse-grid correction technique is used to eliminate the high frequency error caused 
by the data interpolation from the fine-grid subdomain to the coarse-grid subdomain. 

Both driven cavity and jet flow demonstrate the versatility of the proposed multi-grid method. 
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Fig. 1 Configuration of domain decomposition 



Fig. 2b Element layout of driven cavity flow 
for Re = 400 
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